Half-integer Mott-insulator phases in the imbalanced honeycomb lattice 
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Using mean-field theory, we investigate the ground state properties of ultracold bosons loaded in 
a honeycomb lattice with on-site repulsive interactions and imbalanced nearest-neighbor hopping 
amplitudes. Taking into account correlations between strongly coupled neighboring sites through 
an improved Gutzwiller ansatz, we predict the existence of half-integer Mott-insulator phases, i.e. 
states with half-integer filling and vanishing compressibility. These new insulating phases result 
from the interplay between quantum correlations and the topology of the honeycomb lattice, and 
could be easily addressed experimentally as they have clear signatures in momentum space. 



Since its first experimental observation, graphene and 
its remarkable low-energy excitations have sparked a 
whole new branch of research comprising massless Dirac 
fermions, the anomalous quantum Hall effect and spin 
liquids to cite a few [lH3| . These discoveries have found a 
large echo in the field of ultracold atoms loaded in optical 
potentials 0,11, be it for bosons 0,0 or fermions 18, 91. 
In this Letter, we address the bosonic Mott-insulator to 
superfluid (MI-SF) transition in the honeycomb lattice. 
Though widely studied for the square lattice fTol-[l5| , we 
argue here that this transition is richer for graphene be- 



cause of its peculiar topology 16]. Indeed, contrary to 
the square lattice, the honeycomb lattice is a genuine bi- 
partite lattice made of two shifted triangular sub-lattices, 
one labeled with A sites and the other with B sites. The 
honeycomb Bravais unit cell thus possesses two sites and 
one can define symmetric and antisymmetric states in 
the unit cell. This has dramatic consequences for the 
ground state of the interacting system, in either the Mott 
or the superfluid phases. Strikingly, the repulsive Bose- 
Hubbard model on the honeycomb lattice features Mott 
phases at half-integer fillings when nearest-neighbor hop- 
ping amplitudes are imbalanced. In the following, we will 
present our model Hamiltonian and explain how to im- 
prove the usual mean- field approach 17-19) to correctly 
capture the crucial inter-sublattice correlations responsi- 
ble for the half-integer Mott lobes. We will next compute 
the critical value of the MI-SF transition as a function 
of the imbalanced hopping parameter. Finally, we will 
show that these half-filled Mott phases have clear signa- 
tures in momentum space, making them easy to identify 
experimentally. 

In the graphene lattice, each A site is connected to its 
three adjacent B sites by three vectors c a (a = 1,2,3) 
summing up to zero [H| . Their common length define the 
lattice constant a that we set here to unity for conve- 
nience. We investigate the situation where the hopping 
parameter J' along ci -links can be tuned and made dif- 



ferent from the two other identical hopping parameters 
J along C2 and C3 -links. The limit J' « J corresponds 
to weakly-coupled ID-chains, whereas the limit J' ^> J 
corresponds to weakly-coupled dimers. It is worth notic- 
ing that hopping imbalance in the graphene lattice has 
already been achieved experimentally [(Ml . Under these 
conditions, the tight-binding Bose-Hubbard Hamiltonian 
[2o| with on-site repulsive interactions reads: 

H = - J' Y^[ a £ b * + b i a *] ~ J \- a t bi ' + b t fCL ^ 

£ (£,£') 

+ j £[< « - 1) + hi (hi - 1)] - „ xx + hi}. 

(i) 

Here a\ (a a) and h\ (be) represent the creation (annihi- 
lation) operators associated with the A and B sites on a 
Ci-link hereafter labelled £. The corresponding number 
operators are h £ a = a\a^ and = b\be. The (positive) 
interaction strength is U and \i is the chemical potential. 
The summations run over all links £ and, in the kinetic 
term, over their nearest-neighbor ci -links £' such that the 
A-site on £ and the B-site on £' are nearest neighbors [H[ • 
For a given ci-link, there are four such neighboring ci- 
links. 

In the following, we investigate the zero-temperature 
phase dia gra m of Eq.(pQ) within a mean- field ap- 
proach [17|, LL8|, |21[ . We mainly restrict our analysis to the 
dimer regime J' > J. As long as J' < 2 J, the band struc- 
ture of the non-interacting case (U = 0) depicts the cel- 
ebrated conical intersections at the Dirac points around 
E = and the system is a semi- metal. At J' = 2 J, 
the two Dirac points merge and the band structure un- 
dergoes a topological metal-insulator transition 0, |22| : 
when J' > 2 J, the band structure consists of two bands 
separated by 2 (J 7 — 2 J). When J< this is simply the 
energy difference obtained between the symmetric and 
antisymmetric dimer states |±) = (|a) ± |b)) /y/2 (en- 
ergy T^O built on each ci-link. These dimer states then 
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give rise to the two preceding bands, each with a width 
4 J independent of J'. In this weak interlink coupling 
regime (or strong dimer regime), we expect the physics 
to be driven by the lower band and the MI-SF phase tran- 
sition to be controlled by the ratio J/U. In a mean- field 
approach, the Mott ground state is then well approxi- 
mated by Yli where |n, £+) is the Fock state with 
n bosons in the symmetric state of link £. This ground 
state cannot be written as a product of on-site states. 
As such, it is beyond reach of the standard Gutzwiller's 
ansatz. This salient feature directly arises from the 2- 
point topology of the graphene lattice and cannot happen 
with the square lattice where a strong imbalance of one 
hopping parameter leads to weakly-coupled ID chains. 
A simple way to improve on the Gutzwiller's ansatz is to 
take into account correlations between neighboring lat- 
tice sites. This can be done in two ways by writing the 
(normalized) ground state as \GS) = Yle K) : 

\£)=J2fi%\^^ m ^)i (with^|/^| 2 = l) (2) 

n,m n,m 

= E4?>. +;*>->/ (with]T| 5 W| 2 = i), (3) 

p,q p,q 

where \n, A;m, b)^ is the Fock state on the link £ with 
n atoms occupying the A site and m atoms the B site 
while |p, +; g, — } £ is the Fock state on the same link £ 
with p atoms occupying the symmetric state |+) and q 
atoms the anti-symmetric state | — ). If both preceding 
equations represent the most general two-site state, and 
therefore fully describes the system on each ci-link, the 
parametrization Eq.(j3|) proves more useful in the limit 
J' ^> J. It is straightforward to compute the transfor- 
mation between the two set of amplitudes fn}m and g£\ 
by noting that the annihilation and creation operators 
for the |±) states are simply d^ ± = (a^ ± b^)/V2. 

The ground state of Eq.(pQ) has been obtained by 
imaginary-time evolution of an initial state with ran- 
dom values of amplitudes and periodic boundary con- 
ditions. As the structure of the corresponding nonlinear 
time-dependent equations [Io|, is rather involved, a 
4th-order Runge-Kutta method was necessary. We have 
carefully checked that our results were independent of 
the values of numerical parameters such as the number 
of sites, the time- step, the initial value of the amplitudes 
and the maximum occupancy number in the Fock states. 

Uncoupled dimer s. We consider first the situation 
J = and solve for the on-link dimer problem. Figure [T] 
shows the on-site (left) and symmetric (right) average 
occupation numbers and their variances as a function of 
the chemical potential \i and hopping amplitude J' in 
units of the interaction strength U . For small J'/U val- 
ues, the on-site density depicts the usual Mott plateaus 
at integer fillings. But, when J'/U increases, one clearly 
observes the appearance of new plateaus at half-integer 
fillings. At the same time, plateaus at odd integer values 



2p + 1 start to appear in the symmetric state density. For 
larger J'/U values, each of these plateaus splits in two 
new ones with 2p + 1 and 2p + 2 fillings. Correspond- 
ingly (but not shown here), the occupation number of 
the anti-symmetric state decreases from the same value 
2p + 1 at J' = 0, to an almost vanishing value. The evo- 
lution of the on-site and symmetric variances emphasizes 
that the ground state of the system evolves from a on- 
site Fock state to an (almost) Fock state in the |±) basis 
when J'/U is cranked up. Indeed, for J' « 0, the on-site 
variance is almost zero, whereas the variance in the sym- 
metric state is the largest. When J'/U is increased, the 
on-site variance gets larger, whereas the symmetric one 
almost vanishes. More precisely, for J' « 0, the ground 
state is well approximated by Fock states \n, A;n, b)^, 
whereas for larger value of J', it is well approximated by 
Fock states |p, +; 0, — ) £ . One may however note that, for 
J' = 0.6, the symmetric variance, albeit extremely small, 
is not strictly zero and even increases with the chemical 
potential \i. This points at a slight contamination of the 
ground state at finite J' by the anti-symmetric states due 
to the interacting part of the Hamiltonian (pQ). Indeed, in 
the |=b) basis, only the latter is not diagonal and a careful 
inspection shows that it couples states |p, +;<?, — } £ with 
the same total number of atoms p + and only changes 
p and q by steps of two. As a consequence, the actual 
ground state rather reads: 

\GS) t = a \p,+]0-) l + <* q \p-2q,+]2q-) £ . (4) 

q<p/2 

In the large J'/U limit, a is of the order of unity, whereas 
the other coefficients are small and becoming smaller 
with either increasing q or J'/U. Except for p = and 
p = 1, n\_ is always a bit less than an integer and the 
symmetric variance is never strictly zero, even if the dis- 
crepancies goes to when increasing J'/U. This is clearly 
seen in plot (f) of Fig. [1] at J' = 0.6 (dot-dashed line) 
where the variance deviates from zero only for plateaus 
with symmetric filling p > 2 (solid line). For small J'/U 
values, the ground state simply becomes |n, A; n, b)^, such 
that the coefficients a q are now given, up to a normal- 
ization factor, by ( — l) q /(n — q)\q\. They reach a maxi- 
mum for q = n/2, leading to the same number of bosons 
in the |d=) states. This corresponds to a ground state 
with an even- integer symmetric filling p = 2n in Eq. (|4|). 
The case p = 2n + 1 corresponds to an odd number of 
bosons and Eq. (j4]) is not a ground state of the system for 
J' = 0. However, fort large enough J'/U, the additional 
energy cost for this extra boson is compensated for by 
the lowering of the symmetric state energy. It then be- 
comes more favorable, in a given range of fi/U, to have 
the 2n + 1 bosons filling that state. Of course, the full 
transition when J'/U is increased is slightly more com- 
plicated since, at intermediate values, the ground state 
does not simply corresponds to a pure Fock state, either 
in the |±) basis or in the on-site one. 
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Figure 1. (color online). Phase diagram of Hubbard Hamilto- 
nian (pQ) for J = 0. Left: plots of the on-site density n (a) and 
variance Sn (c) as a function of J' /U and ji/U with their color 
codes. Plot (e) is a cut of (c) for J' /U = 0.6. Right: same as 
left but for the symmetric density n+ and variance £n+. For 
large J' /U values, the on-site density depicts plateaus at half- 
integer fillings, corresponding to integer occupation numbers 
of the symmetric state |+). The bottom plots emphasize that 
the first plateau corresponds to the state |1, +; 0, —) £ . 



Coupled dimers. We now consider the (J f /U, fi/U) 
phase diagram in the symmetric state when the hop- 
ping amplitude J is non zero. As depicted by Fig. [2^ 
and Fig. [2)3, the situation is very much similar to the 
usual MI-SF transition. One finds regions with constant 
symmetric occupation numbers n+ and almost vanish- 
ing variance, corresponding to a Mott state with van- 
ishing compressibility These Mott lobes are sur- 
rounded by a superfluid sea where the compressibility is 
finite. For weak values of J/U, the lower Mott plateaus 
are still well visible, whereas the higher plateaus are al- 
most all smoothed out. This is exemplified in Fig. [2k 
for J' /U = 0.6 and J/U = 0.015 where the variance of 
the first two Mott phases are still small. Between the 
plateaus, the variance displays maxima revealing the su- 
perfluid phases. For /i > [/, n + varies smoothly and 
the system is superfluid. For larger value of J/U, the 
half-integer Mott lobes almost entirely disappear except 
at very low values of J'/U, see Fig. [2}d obtained at 
J/U = 0.04. This is exemplified by Fig. Hi at J' jU = 0.6, 
where both n+ and 5n+ vary smoothly as \i is increased. 
The evolution of the probability amplitudes when 
J increases and J' is fixed, is similar to the usual MI- 
SF scenario in the Gutzwiller's ansatz. When J'/U 
is large, the Mott state is approximated by the Fock 



state |po>+;0, — ) £ and g* = S PPo S q o. At the transition, 
the distribution g* starts broadening, but only in the 
symmetric direction q = and the physics takes entirely 
place in the symmetric subspace. For instance, in the 
limit J' ^> J ^> U , the superfluid phase for a density p is 
described by a coherent state in the symmetric space and 
no bosons in the anti-symmetric space, <S> |0)_. At 

intermediate values of J', the Mott state is slightly con- 
taminated by anti-symmetric contributions. However, we 
have checked that the MI-SF transition scenario remains 
the same: g^ still spreads along the symmetric direction 
and keep a structure similar to (|4|). 



(a) (b) 




Figure 2. (color online). Left column: phase diagram of Hub- 
bard Hamiltonian (pQ) for J/U = 0.015. Right column: same 
as the left column but for J/U = 0.04. (a), (b): occupation 
number n+ as a function of J' /U and ji/U and the color code, 
(c), (d) n+ and 5n+ as a function of \i/U for J' /U — 0.6. For 
J — 0.015, the first three Mott plateaus are clearly visible 
whereas they are smoothed out for J — 0.04. 

Critical hopping amplitude. We obtain the MI-SF 
critical hopping rate J C (J' /U, fJ>/U) by monitoring the 
tip of the two first Mott lobes, see Fig. 03 Within 
the usual Bose-Hubbard model, one predicts U / J& = 
z(2p+l + 2a/p(/?+ 1)) where z is the coordination num- 
ber of the lattice [13, [H| and where the integer p is the 
average density. For small J' ' /U the system is made of 
almost independent ID chains (z = 2) and the physics is 
driven by on-site Fock states. The main Mott lobe is at 
p = n = 1 (equivalent ly n + = 2). For imbalanced hop- 
ping parameters, one gets (2 J + J') = [7/5.8 and thus 
J c = J^ — J' ' i ^2. This prediction correctly reproduces 
our numerical results for small values of J' /U . However, 
in the large J' / J limit, the physics takes place in the 
symmetric subspace. It is easy to see from Eq. (pQ) that 
each state |+) is connected to z = 4 nearest-neighbors 
with an effective hopping parameter J/2 giving rise to 
a non-interacting band with finite width 4 J independent 
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of J' . Furthermore, with the help of Eqs (pQ) and ((4]), 
a second-order perturbation shows that, in the large- J 7 
limit, J c = J^P (l — |-jt), in good agreement with our 
numerical data for p = n + = 1 (equivalent ly n = 0.5). 
Remarkably, for large J'/U, the Mott lobe at n = 1 
(n + = 2) saturates, as it should, at the value J c ^ (z = 4, 
p = n + = 2). Finally, even though the present ansatz 
favors one link over the two others, it is quite remark- 
able that our numerical results for the Mott lobe at 
p = n = 1 (n + = 2) cross the critical value J^l\ cor- 
responding to the balanced honeycomb lattice, roughly 
when J c = J' . From that point of view, the system un- 
dergoes a cross-over from a quasi- ID situation with two 
neighbors (weakly-coupled chains) to a 2D situation with 
four neighbors, the balanced honeycomb lattice with 3 
neighbors being the intermediate situation. A thorough 
study of this crossover is beyond the scope of this Letter 
and will be published elsewhere. 
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Figure 3. (color online) Critical value J c for the MI-SF tran- 
sition as a function of J ; (in units of U). The red (resp. 
blue) dots are the numerical values obtained with the ex- 
tended Gutzwiller's ansatz for the n — 1 (resp. n — 0.5) 
Mott lobe. Continuous and long-dashed lines are analytical 
predictions. The horizontal black dashed lines correspond to 
the critical values Jcl\ z = 1, 2, 3. See text for details. 



Experimental signatures. A clear signature of the MI- 
SF transition can be observed in the velocity distribu- 
tion of the atoms when the optical lattice is rapidly 
switched off 



121 |24[. The velocity distribution is then 
simply the momentum distribution in the ground state, 
oc j(GS\c\cj\GS) exp [ik • (Rj — R^)], where R^ is 
the position of the site i. In the usual Guztwiller's ansatz, 
there is no correlation between different sites and only 
terms like (c\ci) or (c\)(cj) exp(zk • (Rj — R^) contribute 
to the signal. With the extended ansatz, we obtain 
additional contributions from terms like (a]^)exp(ik • 
Ci). These nearest-neighbor correlations result in an 
additional modulation of the velocity distribution, a 
smoking-gun of a Mott phase build on a symmetric state. 
For the Fock state |p, +;0, — ), one obtains nM/(k) = 
p (1 + cos (k • ci)). A more accurate description using 



Eq. (j4|) would show that the modulation amplitude is a 
direct measurement of the purity of the ground-state. For 
the superfluid state, a product of on-site coherent states, 
one finds usf0<) = Po^Zij e* k, ( R i -Ri ) (assuming a uni- 
form density po), which depicts peaks at the reciprocal 
lattice vectors modulated by the square of the structure 
factor of the lattice. The additional on-link correlations 
will thus show up as an additional modulation on the top 
of this ideal distribution. The latter properties are con- 
firmed by our numerical calculations, shown in Fig. [4j 
where, choosing ci along Ox, both the superfluid (left 
plot) or the Mott (right plot) phases display a periodic 
modulation along k x (stripes) with period 2tt. Note that, 
beside the preceding modulation, the effect of the struc- 
ture factor in the superfluid phase is clearly visible in the 
different peak heights. 
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Figure 4. (color online) Velocity distribution for J/U = 0.015, 
J'/U = 0.6. Left frame: superfluid phase (/n/U = —0.14). 
Right frame: Mott-insulator phase (/n/U = —0.5). The 
hexagon is the first Brillouin zone with its reciprocal lattice 
vectors (red arrows). Vector ci is chosen along Ox. The 
stripes along k x (with a period 2tt shown by the white arrow) 
are a signature of the quantum correlations between the A 
and B sites along ci -links. See text for details. 

In conclusion, using an extended Gutzwiller's ansatz, 
we have described the properties of the MI-SF transi- 
tion of ultracold bosons in a honeycomb lattice. We 
have found Mott phases at half-integer fillings, arising 
directly from the interplay between quantum correlations 
and the topology of the honeycomb lattice. Future work 
will address the excitations of the system as they can 
lead to additional experimental signatures [25]. Finally, 
it would be interesting to study the impact of an external 
(non-Abelian) gauge field on the properties of the ground 
state [26l-[28| . The Centre for Quantum Technologies is a 
Research Centre of Excellence funded by the Ministry of 
Education and National Research Foundation of Singa- 
pore. This work has been supported by the CNRS-CQT 
LIA FSQL. ChM is a Fellow of the Institute of Advanced 
Studies (NTU). 
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